R Analyses and Figures for Manuscript
Zachary Rowson
6/7/2022
Introduction
This document will contain all analyses and figures produced in R to be included in the manuscript. Note that these are not all the figures to be included, some were produced in PowerPoint and have been saved as image files.
# Raw data, endpoint data, and row objects for tcplfit2
load("../Data/Padilla_DNT60_lmr0_w_egid.Rdata")
load("~/Desktop/StephanieProjects/DNT60Analysis/pipelining/pipelined data/Padilla_DNT60_mc0.rda")
load("../Data/Padilla_DNT60_mc0_n.rda")
load("../Data/Padilla_DNT60_rows_n.rda")
# Results of curve fitting
load("../Data/Padilla_DNT60_tcplfits.rda")
load("../Data/Padilla_DNT60_tcpl_out.Rdata")
tcpl_out.dt <- as.data.table(DNT60_tcpl_out)
# Summarizing results of curve-fitting
load("../Data/Padilla_DNT60_Strong Effectors by Endpoint.Rdata")
load("../Data/Padilla_DNT60_BMC ranges.Rdata") # This data needs to be recompiled
# Results of literature review
load("../Data/Padilla_DNT60_Behavioral LOELs.Rdata")
load("../Data/Literature Review Summary.Rdata")
load("../Data/Literature Review Concentrations and Activity.Rdata")
load("../Data/Literature Review LOELs and Dev LOELs.Rdata") Example Locomotor Response Data
lmr0.egid[cpid=="Chlorpyrifos (ethyl)", cpid := "Chlorpyrifos"]Example tcpl Curve-Fit
Results
Appearance of Endpoint Data
Figure 3: Raw Fluoxetine Vehicle Control Speed Data with Average Speed Endpoint Data Raw and Transformed
# Time Series for Fluoxetine vehicle control
# Units
unit.t = "2 min"
unit.mov = "cm"
unit.conc = paste0("\U03BC","M")
no.A = 10
# Identify movement columns of interest and extract
t.cols <- grep("vt", names(lmr0.egid), value = TRUE)
cols <- t.cols[(no.A+1):length(t.cols)]
A.cols <- t.cols[!(t.cols%in%cols)]
to.fit <- lmr0.egid[wllt=="v" & egid=="E6", -A.cols, with=FALSE]
# Melt data and designate Light and Dark time periods
to.fit_melt <- melt(to.fit[,-"egid"], id.vars = names(to.fit)[1:9], variable.name = "t", value.name = "speed")
to.fit_melt[, t := as.factor( gsub("vt","",t) )]
to.fit_melt[t%in%11:30, phase := "Light"]
to.fit_melt[t%in%31:50, phase := "Dark"]
# Create title, x- and y-axis titles, legend labels, and legend title
title.t <- paste0("Time (",unit.t,")")
label.y <- "Speed"
title.mean <- paste0(label.y, " (",unit.mov,"/",unit.t,")")
title.legend <- "Experimental Phase"
## Create x-axis breaks and labels
m <- 50
x.breaks <- c(11, 20, 30, 40, 50)
x.labels1 <- 2*seq(from=no.A,to=m,by=10)
x.labels2 <- x.labels1 - 2
x.labels <- paste(x.labels2, x.labels1, sep="-")
# Printing of graphic isn't working very well so insert new breaks to make space for edges
x.breaks <- c(5, x.breaks, 55)
x.labels <- c("a", x.labels, "b")
## plot
FlxTSeriesPlot <- ggplot(to.fit_melt) +
geom_boxplot(aes(x=t,y=speed,fill=phase)) +
scale_x_discrete(breaks = x.breaks, labels = x.labels) +
scale_fill_manual(values=c("darkgrey","white")) +
labs(x = title.t, y = title.mean, fill = "Experimental Phase") +
theme_bw() +
theme(text = element_text(size = 26))# Raw endpoint data for Fluoxetine vehicle control
sample_endpoints <- lapply(mc0, function(table) {
table.egid <- gabi::data_egids(table)
table.egid[egid=="E6"]
})
sample_endpoints[c("AUC_L","AUC_D","AUC_T")] = NULL
for (i in 1:13) {
sample_endpoints[[i]][, endp := names(sample_endpoints)[i]]
}
endp.data <- do.call('rbind', sample_endpoints)
max.conc <- max( endp.data[cpid=="Fluoxetine", conc] )
raw_stats <- endp.data[wllt=="v" | (conc==max.conc & cpid=="Fluoxetine")][endp == "avgS_L"][, .(mean=mean(rval,na.rm=T),median=median(rval,na.rm=T)), by=.(conc)]
raw_stats_l <- data.table::melt(raw_stats, id.vars = "conc", variable.name = "stat")
FlxAvgS_LPlot <- ggplot(endp.data[wllt=="v" | (conc==max.conc & cpid=="Fluoxetine")][endp == "avgS_L"]) +
geom_density(aes(x=rval, fill=as.factor(conc), color = as.factor(conc)), alpha = 0.5) +
theme_bw() +
theme(text = element_text(size = 26)) +
scale_color_manual(values = colors[c(1,5)]) +
scale_fill_manual(values = colors[c(1,5)]) +
labs(color = paste0("Concentration (", unit.conc, ")"),
fill = paste0("Concentration (", unit.conc, ")"),
linetype = "",
x = "Raw Endpoint Value",
y = "Density")
FlxAvgS_LPlot <- FlxAvgS_LPlot +
geom_vline(data = raw_stats_l, aes(xintercept=value, linetype=stat), color = colors[c(1,5,1,5)], size=1) +
scale_linetype_discrete(labels = c("Mean","Median")) +
scale_y_continuous(labels = c("0.0","0.1","0.2","0.3","")) +
labs(color = paste0("Concentration (", unit.conc, ")"),
fill = paste0("Concentration (", unit.conc, ")"),
linetype = "",
x = "Raw Endpoint Value",
y = NULL)
legend <- cowplot::get_legend(FlxAvgS_LPlot)
FlxAvgS_LPlot <- FlxAvgS_LPlot + theme(legend.position="none")sample_endpoints_n <- lapply(mc0_n, function(list) {
table <- list[[1]]
table.egid <- gabi::data_egids(table)
table.egid[egid=="E6"]
})
sample_endpoints_n[c("AUC_L","AUC_D","AUC_T")] = NULL
for (i in 1:13) {
sample_endpoints_n[[i]][, endp := names(sample_endpoints_n)[i]]
}
endp.data_n <- do.call('rbind', sample_endpoints_n)
trans_stats <- endp.data_n[wllt=="v" | (conc==max.conc & cpid=="Fluoxetine")][endp == "avgS_L"][, .(mean=mean(rval,na.rm=T),median=median(rval,na.rm=T)), by=.(conc)]
trans_stats_l <- data.table::melt(trans_stats, id.vars = "conc", variable.name = "stat")
FlxAvgS_LPlot_n <-ggplot(endp.data_n[wllt=="v" | (conc==max.conc & cpid=="Fluoxetine")][endp == "avgS_L"]) +
geom_density(aes(x=rval, fill=as.factor(conc), color = as.factor(conc)), alpha = 0.5) +
theme_bw() +
theme(text = element_text(size = 26)) +
scale_color_manual(values = colors[c(1,5)]) +
scale_fill_manual(values = colors[c(1,5)]) +
labs(color = paste0("Concentration (", unit.conc, ")"),
fill = paste0("Concentration (", unit.conc, ")"),
linetype = "",
x = "Transformed Endpoint Value",
y = NULL)
FlxAvgS_LPlot_n <- FlxAvgS_LPlot_n +
geom_vline(data = trans_stats_l, aes(xintercept=value, linetype=stat), color = colors[c(1,5,1,5)], size=1) +
scale_linetype_discrete(labels = c("Mean","Median")) +
labs(color = paste0("Concentration (", unit.conc, ")"),
fill = paste0("Concentration (", unit.conc, ")"),
linetype = "",
x = "Transformed Endpoint Value",
y = NULL)
FlxAvgS_LPlot_n <- FlxAvgS_LPlot_n + theme(legend.position="none")FlxEndpPlots <- cowplot::plot_grid(FlxAvgS_LPlot, FlxAvgS_LPlot_n,
labels = list("B","C"), label_size = 24)
FlxEndpPlots1 <- cowplot::plot_grid(FlxEndpPlots, legend, rel_widths = c(3, .4))
FlxPlots <- cowplot:::plot_grid(FlxTSeriesPlot, FlxEndpPlots1,
labels = list("A"), label_size = 24,
ncol = 1)
FlxPlots
tiff(filename = "test.tiff", height = 5, width = 5, units = "in", res=300)
FlxPlots
dev.off()## png
## 2Supplemental Figure 1: Fluoxetine Vehicle Control and High Concentration Group Raw Endpoint Data
# Raw endpoint data for Fluoxetine vehicle control
raw_stats1 <- endp.data[wllt=="v" | (conc==max.conc & cpid=="Fluoxetine")][, .(mean=mean(rval,na.rm=T),median=median(rval,na.rm=T)), by=.(endp,conc)]
raw_stats1_l <- data.table::melt(raw_stats1, id.vars = c("endp","conc"), variable.name = "stat")
raw_plot_facet <- ggplot(endp.data[wllt=="v" | (conc==max.conc & cpid=="Fluoxetine")]) +
geom_density(aes(x=rval, fill=as.factor(conc), color = as.factor(conc)), alpha = 0.5) +
theme_bw() +
theme(text = element_text(size = 16)) +
facet_wrap(~endp, scales = "free") +
scale_color_manual(values = colors[c(1,5)]) +
scale_fill_manual(values = colors[c(1,5)])
raw_plot_facet <- raw_plot_facet +
geom_vline(data = raw_stats1_l, aes(xintercept=value, linetype=stat), color = rep(colors[c(1,5)],26), size=1) +
scale_linetype_discrete(labels = c("Mean","Median")) +
labs(title = "Raw Endpoint Data for Fluoxetine Vehicle Control and High Concentration Group",
color = paste0("Concentration (", unit.conc, ")"),
fill = paste0("Concentration (", unit.conc, ")"),
x = "Raw Endpoint Value",
y = "Density")
raw_plot_facet
Supplemental Table 1: Calculate means and variance for endpoints and stick into table.
sample.endp.stats <- endp.data[wllt=="v", .(mean=mean(rval,na.rm=T), variance=stats::var(rval,na.rm=T)), by=.(endp)]
sample.endp.stats[, `:=` (mean = signif(mean,digits=3), variance = signif(variance,digits=3))]
DT::datatable(sample.endp.stats, colnames=c("Endpoint Abbreviation","Mean","Variance"),
caption = "Raw Endpoint Data Statistics for Fluoxetine Vehicle Control")Supplemental Table 2: Box-Cox power transformation parameters.
bxcx.params <- lapply(mc0_n, function(data) {
as.data.table( data[c(2,3)] )
})
bxcx.params.dt <- do.call('rbind', bxcx.params)
bxcx.params.dt[, endp := names(bxcx.params)]
datatable(bxcx.params.dt[, .(endp,lam.hat,shift)], colnames=c("Endpoint Abbreviation","Lambda","Shift"),
caption = "Box-Cox Power Transformation Parameters for Endpoints")Supplemental Figure 2: Fluoxetine Vehicle Control and High Concentration Group Transformed Endpoint Data
trans_stats1 <- endp.data_n[wllt=="v" | (conc==max.conc & cpid=="Fluoxetine")][, .(mean=mean(rval,na.rm=T),median=median(rval,na.rm=T)), by=.(endp,conc)]
trans_stats1_l <- data.table::melt(trans_stats1, id.vars = c("endp","conc"), variable.name = "stat")
trans_plot_facet <- ggplot(endp.data_n[wllt=="v" | (conc==max.conc & cpid=="Fluoxetine")]) +
geom_density(aes(x=rval, fill=as.factor(conc), color = as.factor(conc)), alpha = 0.5) +
theme_bw() +
theme(text = element_text(size = 16)) +
facet_wrap(~endp, scales = "free") +
scale_color_manual(values = colors[c(1,5)]) +
scale_fill_manual(values = colors[c(1,5)])
trans_plot_facet <- trans_plot_facet +
geom_vline(data = trans_stats1_l, aes(xintercept=value, linetype=stat), color = rep(colors[c(1,5)],26), size=1) +
scale_linetype_discrete(labels = c("Mean","Median")) +
labs(title = "Raw Endpoint Data for Fluoxetine Vehicle Control and High Concentration Group",
color = paste0("Concentration (", unit.conc, ")"),
fill = paste0("Concentration (", unit.conc, ")"),
x = "Raw Endpoint Value",
y = "Density")
trans_plot_facet
Applying tcplFit2
How many chemical-endpoint pair were active? How many chemicals had at least one active endpoint?
# Endpoints to exclude
exclude <- c("AUC_L","AUC_D","AUC_T")
# Number of active chemical-endpoint pairs
tcpl_out.dt[hitcall>0.8 & !endp%in%exclude, .N]## [1] 57# Number of active chemicals
length( tcpl_out.dt[hitcall>0.8 & !endp%in%exclude, unique(name)] )## [1] 22Evaluating popular fits.
# Popular fits for all endpoint data
tcpl_out.dt[!endp%in%exclude, .N, by = .(fit_method)]## fit_method N
## 1: exp4 272
## 2: poly1 412
## 3: exp5 16
## 4: hill 5
## 5: gnls 47
## 6: pow 38
## 7: exp2 3# Popular fits for active curve-fits
tcpl_out.dt[!endp%in%exclude & hitcall>0.8, .N, by = .(fit_method)]## fit_method N
## 1: exp4 11
## 2: poly1 22
## 3: exp5 5
## 4: gnls 6
## 5: pow 10
## 6: exp2 1
## 7: hill 2Evaluate the the occurrence of developmental toxicity for chemicals fit optimally by the poly1 function.
# Isolate chemicals with active endpoints optimally fitted by poly1
chemicals <- tcpl_out.dt[hitcall>0.8 & fit_method=="poly1" & !endp%in%c("AUC_L","AUC_D","AUC_T"), unique(name)]
# Is developmental toxicity associated with these chemicals?
litReviewConc[, dev_LOEL := min(.SD[devtox==1,conc]), by=.(cpid)]
litReviewConc[is.infinite(dev_LOEL), dev_LOEL := NA]
litReviewConc[cpid%in%chemicals, .(dev_LOEL = unique(dev_LOEL)), by = .(cpid)]## cpid dev_LOEL
## 1: 5-Fluorouracil NA
## 2: Amphetamine 120.0
## 3: Bisphenol A (BPA) 120.0
## 4: Chlorpyrifos 12.0
## 5: Dieldrin 0.4
## 6: Diethylstilbesterol 10.0
## 7: D-sorbitol NA
## 8: Fluoxetine 12.0
## 9: Heptachlor 4.0
## 10: Loperamide 120.0
## 11: Polybrominated diphenyl ether (PBDE)-47 13.0
## 12: Tebuconazole 40.0Activity in Terms of Endpoints
Table 2: Number of chemicals active in each endpoint.
count.hits.endp <- tcpl_out.dt[hitcall>0.8 & !endp%in%exclude, unique(name), by = .(endp)][, .N, by = .(endp)]
datatable(count.hits.endp[order(-N)], colnames=c("Endpoint Abbreviation","Number of Active Chemicals"), caption = "Number of Chemicals Active in Each Endpoint",
rownames = FALSE)Figure 4: Number of chemicals active for each endpoint category and experimental phase.
tcpl_out.actives <- tcpl_out.dt[hitcall>0.8 & !endp%in%c("AUC_L","AUC_D","AUC_T","AUC_r")]
tcpl_out.actives[endp%in%c("avgS_L","avgS_D","avgS_T"), cat1 := "Average Speed"]
tcpl_out.actives[endp%in%c("hbt1_L","hbt1_D","hbt2_L","hbt2_D"), cat1 := "Habituation"]
tcpl_out.actives[endp%in%c("RoA_L","RoA_D"), cat1 := "Range of Activity"]
tcpl_out.actives[endp%in%c("strtlA","strtlAavg","strtlF"), cat1 := "Startle Response"]
tcpl_out.actives[endp%in%c("avgS_L","hbt1_L","hbt2_L","RoA_L"), cat2 := "Light"]
tcpl_out.actives[endp%in%c("avgS_D","hbt1_D","hbt2_D","RoA_D"), cat2 := "Dark"]
tcpl_out.actives[endp%in%c("strtlA","strtlAavg","strtlF"), cat2 := "Transition"]
tcpl_out.actives[endp%in%c("avgS_T"), cat2 := "Total"]
actv.chem.cat1 <- tcpl_out.actives[, .(cat1 = factor(unique(cat1),levels=c("Average Speed","Habituation",
"Range of Activity","Startle Response"))),
by=.(name)]
actv.chem.cat2 <- tcpl_out.actives[, .(cat2 = factor(unique(cat2),levels=c("Light","Transition","Dark","Total"))),
by=.(name)]
activesCatCntPlot <- ggplot(actv.chem.cat1, aes(x=cat1, color=cat1, fill=cat1)) +
geom_bar() +
stat_count(geom = "text", colour = "black", size = 3.5,
aes(label = paste(..count..,"22",sep="/")), position=position_stack(vjust=1.05)) +
theme_bw() +
theme(text = element_text(size=14)) +
labs(x="Endpoint Categories",
y="Number of Chemicals Active in Category",
color="Endpoint Category",
fill="Endpoint Category") +
scale_y_continuous(breaks=c(2,4,6,8,10)) +
scale_x_discrete(labels=c("Average Speed","Habituation","Range of Activity","Startle Response")) +
scale_color_manual(values=viridis::viridis(4),
labels=c("Average Speed","Habituation","Range of Activity","Startle Response")) +
scale_fill_manual(values=viridis::viridis(4),
labels=c("Average Speed","Habituation","Range of Activity","Startle Response"))
activesPhsCntPlot <- ggplot(actv.chem.cat2, aes(x=cat2, color=cat2, fill=cat2)) +
geom_bar() +
stat_count(geom = "text", colour = "black", size = 3.5,
aes(label = paste(..count..,"22",sep="/")), position=position_stack(vjust=1.05)) +
theme_bw() +
theme(text = element_text(size=14)) +
labs(x="Experimental Phase",
y="Number of Chemicals Active in Phase",
color="Experimental Phase",
fill="Experimental Phase") +
scale_y_continuous(breaks=c(2,4,6,8,10,12)) +
scale_color_manual(values=viridis::viridis(4)) +
scale_fill_manual(values=viridis::viridis(4))
cowplot:::plot_grid(activesCatCntPlot, activesPhsCntPlot,
labels = "AUTO", label_size = 24,
ncol = 2)
Supplemtal Figure 3: Combinations of Endpoint Category and Experimental Phase Activity
endp.lists <- split(actv.chem.cat1[,cat1], actv.chem.cat1[,name])
to.plot <- data.table(name = names(endp.lists), list = endp.lists)
phase.lists <- split(actv.chem.cat2[,cat2], actv.chem.cat2[,name])
to.plot1 <- data.table(name = names(phase.lists), list = phase.lists)
plot <- ggplot(to.plot, aes(x=list)) +
geom_bar() +
theme_bw() +
theme(text = element_text(size=14)) +
scale_x_upset(reverse = TRUE) +
labs(x = "Endpoint Category Combinations",
y = "Number of Chemicals with Combination",
title = "Combinations of Endpoint Activity by Endpoint Category")
plot1 <- ggplot(to.plot1, aes(x=list)) +
geom_bar() +
theme_bw() +
theme(text = element_text(size=14)) +
ggupset::scale_x_upset(reverse = TRUE) +
labs(x = "Experimental Phase Combinations",
y = "Number of Chemicals with Combination",
title = "Combinations of Endpoint Activity by Experimental Phase")
cowplot::plot_grid(plot, plot1, rel_widths = c(9,8),
labels = "AUTO", label_size = 24)
Directionality of endpoint activity was taken from BMC heatmap which will be presented in BMC section
Results in Terms of Chemicals
Chemicals that were active in endpoints other than Average Speed endpoints.
active.avgS <- tcpl_out.dt[hitcall>0.8 & endp%in%c("AUC_L","AUC_D","AUC_T"), unique(name)]
tcpl_out.dt[hitcall>0.8 & !name%in%active.avgS & !endp%in%exclude, unique(name)]## [1] "5,5-Diphenylhydantoin" "5-Fluorouracil" "Acrylamide"
## [4] "Bis(tributyltin) oxide" "Chlorpyrifos" "Diethylstilbesterol"
## [7] "Dieldrin" "D-sorbitol" "Heptachlor epoxide"
## [10] "Phenobarbital" "Tebuconazole" "Triethyltin"Linear Regression of Chlorpyrifos Average Speed in Dark
cpf_row <- rows_n[["Chlorpyrifos"]][["avgS_D"]]
cpf_dt <- data.frame(conc = c(cpf_row$conc,rep(0,length(cpf_row$bresp))),
resp = c(cpf_row$resp,cpf_row$bresp))
cpf_dt$conc = cpf_dt$conc - mean(cpf_dt$conc)
# Regress responses with a linear model
cpf_lm <- lm(resp ~ conc, data = cpf_dt)
confint(cpf_lm)## 2.5 % 97.5 %
## (Intercept) -0.02297693 0.04307083
## conc -0.04239414 0.01099271Figure 5: Behavior profiles of Chlorpyrifos with active endpoints.
chemical <- "Chlorpyrifos"
# Units
unit.t = "min"
unit.mov = "cm"
unit.conc = paste0("\U03BC","M")
prsp = "SA"
no.A = 10
# Extract chemical data
group <- unique(lmr0.egid[cpid == chemical, egid])
## Identify movement columns of interest
t.cols <- grep("vt", names(lmr0.egid), value = TRUE)
cols <- t.cols[(no.A+1):length(t.cols)]
A.cols <- t.cols[!(t.cols%in%cols)]
## extract data to be plotted, exclude acclimation
to.fit <- lmr0.egid[cpid==chemical | (wllt=="v" & egid==group), -A.cols, with=FALSE]
# create appropriate axes titles for plots
label.y <- "Speed"
# Format data for plotting
## calculate mean and 50% CIs for each vector column by concentration group, excluding concentration
exclude.A <- t.cols[!(t.cols%in%A.cols)]
means <- to.fit[, lapply(.SD, function(col) mean(col,na.rm=T)),
.SDcols = exclude.A,
by = conc]
## calculate CI's for transformed values then transform back
shift <- 1
logCIs <- to.fit[, lapply(.SD, function(x) log10(x+shift)), .SDcols=exclude.A, by=conc][
, lapply(.SD, function(x) t.test(x,conf.level=0.50)$conf.int), .SDcols=exclude.A, by=conc]
CIs <- logCIs[, lapply(.SD, function(x) (10^x)-shift), by=conc][,lapply(.SD, function(col) abs(diff(col))/2), .SDcols=exclude.A, by=conc]
## elongate means and CIs data, and join
means_long <- data.table::melt(means, id.vars = "conc", variable.name = "t", value.name = "mean")
means_long[, t := sub("vt","",t)]
CIs_long <- data.table::melt(CIs, id.vars = "conc", variable.name = "t", value.name = "CI")
CIs_long[, t := sub("vt","",t)]
stats <- means_long[CIs_long, on = c("conc","t")][, conc := as.factor(conc)]
stats[, t := as.numeric(t)]
# create standard error of mean estimates by time period and plot as ribbons or error bars
# plot time-series data
## create title, x- and y-axis titles, legend labels, and legend title
title <- paste0("Sample Averaged Time-Series for ", chemical)
title.t <- paste0("Time (",unit.t,")")
title.mean <- paste0("Mean ", label.y, " (",unit.mov,"/",unit.t,")")
conc.n <- to.fit[wllq==1, .N, by=.(conc)][order(conc)]
legend.labels <- paste0(conc.n$conc, ", n=", conc.n$N)
title.legend <- paste0("Concentration (", unit.conc, ")")
## get better colors for plotting
N <- length(unique(to.fit[,conc]))
colors <- viridis::viridis(N)
## create x-axis breaks and labels
m <- as.integer(max(means_long[,t]))
x.breaks <- seq(from=no.A,to=m,by=10)
x.labels1 <- 2*seq(from=no.A,to=m,by=10)
x.labels2 <- x.labels1 - 2
x.labels <- paste(x.labels2, x.labels1, sep="-")
## plot
SAplot_CPF <- ggplot() +
geom_point(data = stats, aes(x=t, y=mean, color=as.factor(conc))) +
geom_line(data = stats, aes(x=t, y=mean, color=conc, group=conc)) +
scale_x_continuous(breaks = x.breaks, labels = x.labels) +
scale_color_manual(values = colors, labels=legend.labels) +
geom_ribbon(data = stats,
aes(x=t, ymax=mean+CI, ymin=mean-CI, group=conc, fill=conc),
alpha = 0.25) +
geom_rect(aes(xmin=10,xmax=30,ymin=-1.25,ymax=-0.5) ,fill="white", color="black") +
geom_rect(aes(xmin=30,xmax=50,ymin=-1.25,ymax=-0.5) ,fill="black", color="black") +
annotate("text", x=c(20,40), y=rep(-0.875,2), color=c("black","white"), label=c("Light Phase","Dark Phase"), size=5) +
scale_fill_manual(values = colors, labels=legend.labels) +
labs(title = title, subtitle = "Acclimation Period Excluded: 50% Confidence Bands",
x = title.t, y = title.mean, color = title.legend) +
guides(fill = "none") +
theme_bw() +
theme(text = element_text(size = 18))row1_CPF <- rows_n[[chemical]][["strtlAavg"]]
row2_CPF <- rows_n[[chemical]][["hbt1_D"]]
fit1_CPF <- concRespCoreZR(row1_CPF, do.plot = TRUE, verbose.plot = FALSE)[["plot"]]
fit2_CPF <- concRespCoreZR(row2_CPF, do.plot = TRUE, verbose.plot = FALSE)[["plot"]]
# fit1_CPF$theme$text$size = 4
# fit2_CPF$theme$text$size = 4
#
# fit1_CPF$heme$line$size = 0.15# grobs <- list(SAplot_CPF, fit1_CPF, fit2_CPF)
# layout <- matrix(c(1,1,1,1,2,3), nrow=3, ncol=2, byrow=TRUE)
#
# gridExtra::grid.arrange(grobs = grobs, layout_matrix = layout)p1 <- cowplot::plot_grid(fit1_CPF, fit2_CPF, nrow=1)
cowplot::plot_grid(SAplot_CPF, p1, nrow=2)
BMC Heatmap to be displayed later.
ac <- c("avgS_L", "hbt1_L", "hbt2_L","RoA_L",
"strtlA","strtlAavg", "strtlF",
"avgS_D","hbt1_D", "hbt2_D", "RoA_D",
"avgS_T","AUC_r")
## Gather BMC's
bmd <- lapply(tcplfits_n, function(chm) unlist(lapply(chm, function(fit) {
hitcall <- fit[["summary"]]$hitcall
bmd <- fit[["summary"]]$bmd
if (!(hitcall>0.8 & !is.na(bmd) & bmd != 0)) {
bmd <- 10000
}
bmd
})))
bmd <- do.call("rbind", bmd)
bmd <- log10(bmd[,colnames(bmd)%in%ac])
# bmd <- bmd[,-c(1,9,14)]
# Identify active chemicals
actives <- rownames(bmd)[apply(bmd, 1, function(row) any(row < 4))]
# Isolate data of interest
to.fit <- bmd
# PBDE-47 name is too long
rownames(bmd)[50] <- "PBDE-47"
actives[19] <- "PBDE-47"# Create a matrix specifying if an up, down, up and down, or down and up arrow should be printed in cells.
to.fit <- lapply(tcplfits_n, function(chm) chm[-c(1,9,14)])
layer.mat <- lapply(to.fit, function(chm) unlist(lapply(chm, function(fit) {
hitcall <- fit[["summary"]]$hitcall
fit_method <- fit[["summary"]]$fit_method
dir <- sign(fit[["summary"]]$top)
if (hitcall>0.8 & fit_method=="gnls" & dir==1) {
layer <- 0
} else if (hitcall>0.8 & fit_method=="gnls" & dir==-1) {
layer <- 1
} else if (hitcall>0.8 & dir==1) {
layer <- 2
} else if (hitcall>0.8 & dir==-1) {
layer <- 3
} else layer <- 4
})))
layer.mat <- do.call("rbind", layer.mat)
rownames(layer.mat)[50] <- "PBDE-47" # Change PBDE-47 name to match BMD matrix
layer.mat1 <- layer.mat[row.names(layer.mat) %in% actives,]# Column labels
col_labels <- c(expression("Average Speed in Light"^1), "Habiutation 1 in Light", "Habituation 2 in Light", expression("Range of Activity in Light"),
"Startle Acceleration", "Startle Relative to Avg. Speed in Light", "Startle Fold-Change",
expression("Average Speed in Dark"^1), "Habituation 1 in Dark", "Habituation 2 in Dark", expression("Range of Activity in Dark"),
"Average Speed in Both Phases", expression("AUC in Dark / AUC in Light Ratio"^2)) # Superscripts notate references in poster
# Legends
# Custom heat legend.
f2 = circlize::colorRamp2(seq(min(bmd), max(bmd), length = 8), rev(viridis(8)), space = "sRGB")
heat_lgd = Legend(col_fun = f2,
title = paste0("BMC log(","\U03BC","M)"),
title_position = "lefttop",
legend_width = unit(4,"cm"),
direction = "horizontal")
# Column annotation by phase legend.
ann_lgd = Legend(labels = c("Light","Transition","Dark","Light+Dark"),
title = "Phase",
title_position = "leftcenter",
labels_gp = gpar(fontsize=8),
title_gp = gpar(fontsize=8),
legend_gp = grid::gpar(fill = c("white","grey","black","red")),
border = TRUE,
nrow = 1,
column_gap = unit(5, 'mm'))
# Legend for cell signal arrows.
dir_lgd = Legend(labels = c(paste("\U2191","Gain"),
paste("\U2193","Loss"),
paste("\U21C5","GainLoss"),
paste("\U21F5","LossGain")),
title = "Signal Direction",
labels_gp = gpar(fontsize=8),
title_gp = gpar(fontsize=8),
title_position = "leftcenter",
nrow = 1,
column_gap = unit(0, 'mm')) # Will produce warnings, don't worry.
lgd_list <- packLegend(ann_lgd, dir_lgd)
# Create column annotation indicating the phase of the LMR that is described.
column_ha <- ComplexHeatmap::HeatmapAnnotation(Phase = factor(c(rep("Light",4), rep("Transition",3),rep("Dark",4),rep("Light+Dark",2)),
levels=c("Light","Transition","Dark","Light+Dark")),
border = TRUE,
simple_anno_size = unit(0.25, 'cm'),
col = list(Phase=c("Light"="white",
"Transition"="grey",
"Dark"="black",
"Light+Dark"="red")),
annotation_legend_param = list(nrow = 1),
show_annotation_name = FALSE,
show_legend = FALSE)
# Create function to add arrows indicating signal direction in cells.
cell_fun <- function(j, i, x, y, width, height, fill) {
if (layer.mat1[i,j] == 0) {
grid.text("\U21C5", x, y, gp=gpar(fontsize=6))
} else if (layer.mat1[i,j] == 1) {
grid.text("\U21F5", x, y, gp=gpar(fontsize=6))
} else if (layer.mat1[i,j] == 2) {
grid.text("\U2191", x, y, gp=gpar(fontsize=6))
} else if (layer.mat1[i,j] == 3) {
grid.text("\U2193", x, y, gp=gpar(fontsize=6))
}
}# Create main heat map.
htlist <- ComplexHeatmap::Heatmap(bmd[row.names(bmd) %in% actives,],
# Specify some parameters for heat legend.
name = paste0("BMC log(","\U03BC","M)"),
col = f2,
border_gp = grid::gpar(col="black",lwd=1),
rect_gp=grid::gpar(col="grey"),
show_heatmap_legend = TRUE,
heatmap_legend_param = list(legend_height = unit(3.5,"cm"),
direction = "vertical",
title_gp = grid::gpar(fontsize=10),
labels_gp = grid::gpar(fontsize=10)),
# Heatmap width and height
# width = unit(16,"in"),
# height = unit(18,"in"),
# Column label parameters
column_labels = col_labels, column_names_rot = 45,
# Split columns by phase.
column_split = factor(c(rep("Light",4), rep("Transition",3),rep("Dark",4),rep("Light+Dark",2)),
levels=c("Light","Transition","Dark","Light+Dark")),
# Specify some column parameters.
column_title = NULL,
top_annotation = column_ha,
# Add signal direction arrows.
cell_fun = cell_fun,
# Specify some parameters row dendrogram aesthetics and row labels.
row_title_side = "right",
row_title_rot = 0,
row_split = 5,
row_dend_side = "right",
row_names_side = "left",
# Clustering parameters.
cluster_columns = FALSE,
cluster_rows = TRUE,
clustering_distance_rows = "pearson",
# Font sizes
row_names_gp = grid::gpar(fontsize=8),
column_names_gp = grid::gpar(fontsize=10),
row_title_gp = grid::gpar(fontsize=10),
column_title_gp = grid::gpar(fontsize=10)
)Supplemental Figure 4: Activity in endpoints other than Average Speed endpoints.
isolate <- actives[order(actives)]
isolate <- c(isolate[c(1,2,6,8,21)], isolate[c(10,12,13,20,22)])
layer.mat2 <- layer.mat[isolate,]
cell_fun1 <- function(j, i, x, y, width, height, fill) {
if (layer.mat2[i,j] == 0) {
grid.text("\U21C5", x, y, gp=gpar(fontsize=12))
} else if (layer.mat2[i,j] == 1) {
grid.text("\U21F5", x, y, gp=gpar(fontsize=12))
} else if (layer.mat2[i,j] == 2) {
grid.text("\U2191", x, y, gp=gpar(fontsize=12))
} else if (layer.mat2[i,j] == 3) {
grid.text("\U2193", x, y, gp=gpar(fontsize=12))
}
}htlist2 <- ComplexHeatmap::Heatmap(bmd[isolate,],
# Specify some parameters for heat legend.
name = paste0("BMC log(","\U03BC","M)"),
col = f2,
border_gp = grid::gpar(col="black",lwd=1),
rect_gp=grid::gpar(col="grey"),
show_heatmap_legend = TRUE,
heatmap_legend_param = list(legend_height = unit(3.5,"cm"),
direction = "vertical",
title_gp = grid::gpar(fontsize=10),
labels_gp = grid::gpar(fontsize=10)),
# Heatmap width and height
# width = unit(16,"in"),
# height = unit(18,"in"),
# Column label parameters
column_labels = col_labels, column_names_rot = 45,
# Split columns by phase and row by activity type
column_split = factor(c(rep("Light",4), rep("Transition",3),rep("Dark",4),rep("Light+Dark",2)),
levels=c("Light","Transition","Dark","Light+Dark")),
# row_split = factor(c(rep("1",5), rep("2",5))),
# Specify some column parameters.
column_title = NULL,
top_annotation = column_ha,
# Add signal direction arrows.
cell_fun = cell_fun1,
# Specify some parameters row dendrogram aesthetics and row labels.
row_title = NULL,
row_title_rot = 0,
#row_split = 3,
row_names_side = "left",
# Clustering parameters.
cluster_columns = FALSE,
cluster_rows = FALSE,
# Font sizes
row_names_gp = grid::gpar(fontsize=8),
column_names_gp = grid::gpar(fontsize=10),
row_title_gp = grid::gpar(fontsize=10),
column_title_gp = grid::gpar(fontsize=10)
)
ht <- draw(htlist2, merge_legend = FALSE, annotation_legend_list = lgd_list,
annotation_legend_side = "top", align_annotation_legend = "heatmap_center")
ro <- row_order(ht)
co <- column_order(ht)
# Highlight startle endpoints and average acceleration in Dark endpoint
decorate_heatmap_body(paste0("BMC log(","\U03BC","M)"), row_slice = 1, column_slice = 1, {
grid.rect(unit(1.025,'npc'), unit(0.7,'npc'),
width = (length(co[[1]]) + length(co[[2]]))/length(co[[1]]) * unit(0.415,'npc') + unit(1,'mm'),
height = (length(ro[[1]]) + length(ro[[2]]))/length(ro[[1]]) * unit(0.344,'npc') + unit(1,'mm'),
gp = gpar(lwd=2, lty=1, fill=NA, col="red"), just=c('left','top')
)
})
decorate_heatmap_body(paste0("BMC log(","\U03BC","M)"), row_slice = 1, column_slice = 1, {
grid.rect(unit(2.05,'npc'), unit(1,'npc'),
width = (length(co[[1]]) + length(co[[2]]))/length(co[[1]]) * unit(0.13,'npc') + unit(1,'mm'),
height = (length(ro[[1]]) + length(ro[[2]]))/length(ro[[1]]) * unit(0.245,'npc') + unit(1,'mm'),
gp = gpar(lwd=2, lty=1, fill=NA, col="red"), just=c('left','top')
)
})
Supplemental Table 3: Number of hits per chemical.
count.hits.chm <- tcpl_out.dt[hitcall>0.8 & !endp%in%exclude, .N, by = .(name)]
datatable(count.hits.chm[order(-N)], colnames=c("Chemical Name","Number of Hits"), caption = "Number of Hits per Chemical",
rownames = FALSE)Supplemental Figure 5: Create plots for comparing Paraquat and Heptachlor epoxide activity profiles.
chemical <- "Heptachlor epoxide"
# Units
unit.t = "min"
unit.mov = "cm"
unit.conc = paste0("\U03BC","M")
prsp = "SA"
no.A = 10
# Extract chemical data
group <- unique(lmr0.egid[cpid == chemical, egid])
## Identify movement columns of interest
t.cols <- grep("vt", names(lmr0.egid), value = TRUE)
cols <- t.cols[(no.A+1):length(t.cols)]
A.cols <- t.cols[!(t.cols%in%cols)]
## extract data to be plotted, exclude acclimation
to.fit <- lmr0.egid[cpid==chemical | (wllt=="v" & egid==group), -A.cols, with=FALSE]
# create appropriate axes titles for plots
label.y <- "Speed"
# Format data for plotting
## calculate mean and 50% CIs for each vector column by concentration group, excluding concentration
exclude.A <- t.cols[!(t.cols%in%A.cols)]
means <- to.fit[, lapply(.SD, function(col) mean(col,na.rm=T)),
.SDcols = exclude.A,
by = conc]
## calculate CI's for transformed values then transform back
shift <- 1
logCIs <- to.fit[, lapply(.SD, function(x) log10(x+shift)), .SDcols=exclude.A, by=conc][
, lapply(.SD, function(x) t.test(x,conf.level=0.50)$conf.int), .SDcols=exclude.A, by=conc]
CIs <- logCIs[, lapply(.SD, function(x) (10^x)-shift), by=conc][,lapply(.SD, function(col) abs(diff(col))/2), .SDcols=exclude.A, by=conc]
## elongate means and CIs data, and join
means_long <- data.table::melt(means, id.vars = "conc", variable.name = "t", value.name = "mean")
means_long[, t := sub("vt","",t)]
CIs_long <- data.table::melt(CIs, id.vars = "conc", variable.name = "t", value.name = "CI")
CIs_long[, t := sub("vt","",t)]
stats <- means_long[CIs_long, on = c("conc","t")][, conc := as.factor(conc)]
stats[, t := as.numeric(t)]
# create standard error of mean estimates by time period and plot as ribbons or error bars
# plot time-series data
## create title, x- and y-axis titles, legend labels, and legend title
title <- paste0("Sample Averaged Time-Series for ", chemical)
title.t <- paste0("Time (",unit.t,")")
title.mean <- paste0("Mean ", label.y, " (",unit.mov,"/",unit.t,")")
conc.n <- to.fit[wllq==1, .N, by=.(conc)][order(conc)]
legend.labels <- paste0(conc.n$conc, ", n=", conc.n$N)
title.legend <- paste0("Concentration (", unit.conc, ")")
## get better colors for plotting
N <- length(unique(to.fit[,conc]))
colors <- viridis::viridis(N)
## create x-axis breaks and labels
m <- as.integer(max(means_long[,t]))
x.breaks <- seq(from=no.A,to=m,by=10)
x.labels1 <- 2*seq(from=no.A,to=m,by=10)
x.labels2 <- x.labels1 - 2
x.labels <- paste(x.labels2, x.labels1, sep="-")
## plot
plot.HepEpox <- ggplot() +
geom_point(data = stats, aes(x=t, y=mean, color=as.factor(conc))) +
geom_line(data = stats, aes(x=t, y=mean, color=conc, group=conc)) +
scale_x_continuous(breaks = x.breaks, labels = x.labels) +
scale_color_manual(values = colors, labels=legend.labels) +
geom_ribbon(data = stats,
aes(x=t, ymax=mean+CI, ymin=mean-CI, group=conc, fill=conc),
alpha = 0.25) +
geom_rect(aes(xmin=10,xmax=30,ymin=-1.25,ymax=-0.5) ,fill="white", color="black") +
geom_rect(aes(xmin=30,xmax=50,ymin=-1.25,ymax=-0.5) ,fill="black", color="black") +
annotate("text", x=c(20,40), y=rep(-0.875,2), color=c("black","white"), label=c("Light Phase","Dark Phase"), size=3) +
scale_fill_manual(values = colors, labels=legend.labels) +
labs(title = title, subtitle = "Acclimation Period Excluded: 50% Confidence Bands",
x = title.t, y = title.mean, color = title.legend) +
guides(fill = "none") +
theme_bw() +
theme(text = element_text(size = 14))signal.HepEpox <- tcpl_out.dt[name=="Heptachlor epoxide" & hitcall>0.8 & !endp%in%exclude, .(endp,top_over_cutoff)]
table.HepEpox <- gridExtra::tableGrob(signal.HepEpox, rows = NULL, cols = c("Endpoint Abbreviation", "Top Over Cutoff"))chemical <- "Paraquat"
# Units
unit.t = "min"
unit.mov = "cm"
unit.conc = paste0("\U03BC","M")
prsp = "SA"
no.A = 10
# Extract chemical data
group <- unique(lmr0.egid[cpid == chemical, egid])
## Identify movement columns of interest
t.cols <- grep("vt", names(lmr0.egid), value = TRUE)
cols <- t.cols[(no.A+1):length(t.cols)]
A.cols <- t.cols[!(t.cols%in%cols)]
## extract data to be plotted, exclude acclimation
to.fit <- lmr0.egid[cpid==chemical | (wllt=="v" & egid==group), -A.cols, with=FALSE]
# create appropriate axes titles for plots
label.y <- "Speed"
# Format data for plotting
## calculate mean and 50% CIs for each vector column by concentration group, excluding concentration
exclude.A <- t.cols[!(t.cols%in%A.cols)]
means <- to.fit[, lapply(.SD, function(col) mean(col,na.rm=T)),
.SDcols = exclude.A,
by = conc]
## calculate CI's for transformed values then transform back
shift <- 1
logCIs <- to.fit[, lapply(.SD, function(x) log10(x+shift)), .SDcols=exclude.A, by=conc][
, lapply(.SD, function(x) t.test(x,conf.level=0.50)$conf.int), .SDcols=exclude.A, by=conc]
CIs <- logCIs[, lapply(.SD, function(x) (10^x)-shift), by=conc][,lapply(.SD, function(col) abs(diff(col))/2), .SDcols=exclude.A, by=conc]
## elongate means and CIs data, and join
means_long <- data.table::melt(means, id.vars = "conc", variable.name = "t", value.name = "mean")
means_long[, t := sub("vt","",t)]
CIs_long <- data.table::melt(CIs, id.vars = "conc", variable.name = "t", value.name = "CI")
CIs_long[, t := sub("vt","",t)]
stats <- means_long[CIs_long, on = c("conc","t")][, conc := as.factor(conc)]
stats[, t := as.numeric(t)]
# create standard error of mean estimates by time period and plot as ribbons or error bars
# plot time-series data
## create title, x- and y-axis titles, legend labels, and legend title
title <- paste0("Sample Averaged Time-Series for ", chemical)
title.t <- paste0("Time (",unit.t,")")
title.mean <- paste0("Mean ", label.y, " (",unit.mov,"/",unit.t,")")
conc.n <- to.fit[wllq==1, .N, by=.(conc)][order(conc)]
legend.labels <- paste0(conc.n$conc, ", n=", conc.n$N)
title.legend <- paste0("Concentration (", unit.conc, ")")
## get better colors for plotting
N <- length(unique(to.fit[,conc]))
colors <- viridis::viridis(N)
## create x-axis breaks and labels
m <- as.integer(max(means_long[,t]))
x.breaks <- seq(from=no.A,to=m,by=10)
x.labels1 <- 2*seq(from=no.A,to=m,by=10)
x.labels2 <- x.labels1 - 2
x.labels <- paste(x.labels2, x.labels1, sep="-")
## plot
plot.Paraquat <- ggplot() +
geom_point(data = stats, aes(x=t, y=mean, color=as.factor(conc))) +
geom_line(data = stats, aes(x=t, y=mean, color=conc, group=conc)) +
scale_x_continuous(breaks = x.breaks, labels = x.labels) +
scale_color_manual(values = colors, labels=legend.labels) +
geom_ribbon(data = stats,
aes(x=t, ymax=mean+CI, ymin=mean-CI, group=conc, fill=conc),
alpha = 0.25) +
geom_rect(aes(xmin=10,xmax=30,ymin=-1.25,ymax=-0.5) ,fill="white", color="black") +
geom_rect(aes(xmin=30,xmax=50,ymin=-1.25,ymax=-0.5) ,fill="black", color="black") +
annotate("text", x=c(20,40), y=rep(-0.875,2), color=c("black","white"), label=c("Light Phase","Dark Phase"), size=3) +
scale_fill_manual(values = colors, labels=legend.labels) +
labs(title = title, subtitle = "Acclimation Period Excluded: 50% Confidence Bands",
x = title.t, y = title.mean, color = title.legend) +
guides(fill = "none") +
theme_bw() +
theme(text = element_text(size = 14))signal.Paraquat <- tcpl_out.dt[name=="Paraquat" & hitcall>0.8 & !endp%in%exclude, .(endp,top_over_cutoff)]
table.Paraquat <- gridExtra::tableGrob(signal.Paraquat, rows = NULL, cols = c("Endpoint Abbreviation", "Top Over Cutoff"))cowplot::plot_grid(plot.HepEpox, plot.Paraquat, table.HepEpox, table.Paraquat, ncol = 2, nrow = 2,rel_heights = c(3,1))
Potential Reference Chemicals
Supplemental Table 4: Active Chemicals Inducing the Strongest Signal in Endpoints
aenms <- strgEffectors_byEndp[, unique(aenm)]
# Format data and print as a good table
strgEffectors_formatted <- strgEffectors_byEndp[order(aenm), .(aenm,direction,cpid,conc,rval)]
strgEffectors_formatted[, rval:=round(rval, digits=2)]
setnames(strgEffectors_formatted, names(strgEffectors_formatted),
c("Assay Endpoint Name","Signal Direction","Chemical Name",paste0("Concentration ","(\U03BC","M)"),"SEs from Control Mean"))
datatable(strgEffectors_formatted)BMC Estimates
Table 4: Which chemicals appear most potent by median BMC across active endpoints?
select <- apply(bmd, 1, function(row) any(row<4))
med_BMC <- bmd[select,]
med_BMC[med_BMC==4] <- NA
med_BMC.1 <- apply(med_BMC, 1, median, na.rm=TRUE)
med_BMC.2 <- 10^(med_BMC.1[order(med_BMC.1)])
medBMC.dt <- data.table("Chemical Name" = names(med_BMC.2), "Median BMC" = signif(med_BMC.2, digits=3))
datatable(medBMC.dt, rownames = FALSE,
caption = "Median BMCs of Chemicals Across Active Endpoints",
colnames = c("Chemical Name", paste0("Median BMC (","\U03BC","M)")))Table 5: BMC ranges for chemicals.
BMDRanges[, BMCrange := signif(BMDRange, digits=3)]
datatable(BMDRanges, rownames = FALSE,
caption = "BMC Ranges for Chemicals with Multiple Hits",
colnames = c("Chemical Name", paste0("Range of BMCs (","\U03BC","M)")))Supplemental Figure 6: Diazepam behavior data and active curve-fits.
plot_SA <- function(chemical) {
# Units
unit.t = "min"
unit.mov = "cm"
unit.conc = paste0("\U03BC","M")
prsp = "SA"
no.A = 10
# Extract chemical data
group <- unique(lmr0.egid[cpid == chemical, egid])
## Identify movement columns of interest
t.cols <- grep("vt", names(lmr0.egid), value = TRUE)
cols <- t.cols[(no.A+1):length(t.cols)]
A.cols <- t.cols[!(t.cols%in%cols)]
## extract data to be plotted, exclude acclimation
to.fit <- lmr0.egid[cpid==chemical | (wllt=="v" & egid==group), -A.cols, with=FALSE]
# create appropriate axes titles for plots
label.y <- "Speed"
# Format data for plotting
## calculate mean and 50% CIs for each vector column by concentration group, excluding concentration
exclude.A <- t.cols[!(t.cols%in%A.cols)]
means <- to.fit[, lapply(.SD, function(col) mean(col,na.rm=T)),
.SDcols = exclude.A,
by = conc]
## calculate CI's for transformed values then transform back
shift <- 1
logCIs <- to.fit[, lapply(.SD, function(x) log10(x+shift)), .SDcols=exclude.A, by=conc][
, lapply(.SD, function(x) t.test(x,conf.level=0.50)$conf.int), .SDcols=exclude.A, by=conc]
CIs <- logCIs[, lapply(.SD, function(x) (10^x)-shift), by=conc][,lapply(.SD, function(col) abs(diff(col))/2), .SDcols=exclude.A, by=conc]
## elongate means and CIs data, and join
means_long <- data.table::melt(means, id.vars = "conc", variable.name = "t", value.name = "mean")
means_long[, t := sub("vt","",t)]
CIs_long <- data.table::melt(CIs, id.vars = "conc", variable.name = "t", value.name = "CI")
CIs_long[, t := sub("vt","",t)]
stats <- means_long[CIs_long, on = c("conc","t")][, conc := as.factor(conc)]
stats[, t := as.numeric(t)]
# create standard error of mean estimates by time period and plot as ribbons or error bars
# plot time-series data
## create title, x- and y-axis titles, legend labels, and legend title
title <- paste0("Sample Averaged Time-Series for ", chemical)
title.t <- paste0("Time (",unit.t,")")
title.mean <- paste0("Mean ", label.y, " (",unit.mov,"/",unit.t,")")
conc.n <- to.fit[wllq==1, .N, by=.(conc)][order(conc)]
legend.labels <- paste0(conc.n$conc, ", n=", conc.n$N)
title.legend <- paste0("Concentration (", unit.conc, ")")
## get better colors for plotting
N <- length(unique(to.fit[,conc]))
colors <- viridis::viridis(N)
## create x-axis breaks and labels
m <- as.integer(max(means_long[,t]))
x.breaks <- seq(from=no.A,to=m,by=10)
x.labels1 <- 2*seq(from=no.A,to=m,by=10)
x.labels2 <- x.labels1 - 2
x.labels <- paste(x.labels2, x.labels1, sep="-")
## plot
plot <- ggplot() +
geom_point(data = stats, aes(x=t, y=mean, color=as.factor(conc))) +
geom_line(data = stats, aes(x=t, y=mean, color=conc, group=conc)) +
scale_x_continuous(breaks = x.breaks, labels = x.labels) +
scale_color_manual(values = colors, labels=legend.labels) +
geom_ribbon(data = stats,
aes(x=t, ymax=mean+CI, ymin=mean-CI, group=conc, fill=conc),
alpha = 0.25) +
geom_rect(aes(xmin=10,xmax=30,ymin=-1.25,ymax=-0.5) ,fill="white", color="black") +
geom_rect(aes(xmin=30,xmax=50,ymin=-1.25,ymax=-0.5) ,fill="black", color="black") +
annotate("text", x=c(20,40), y=rep(-0.875,2), color=c("black","white"), label=c("Light Phase","Dark Phase"), size=5) +
scale_fill_manual(values = colors, labels=legend.labels) +
labs(title = title, subtitle = "Acclimation Period Excluded: 50% Confidence Bands",
x = title.t, y = title.mean, color = title.legend) +
guides(fill = "none") +
theme_bw() +
theme(text = element_text(size = 18))
return(plot)
}chemical <- "Diazepam"
SAplot_Dzp <- plot_SA(chemical)
row1_Dzp <- rows_n[[chemical]][["avgS_L"]]
row2_Dzp <- rows_n[[chemical]][["RoA_L"]]
row3_Dzp <- rows_n[[chemical]][["avgS_T"]]
fit1_Dzp <- concRespCoreZR(row1_Dzp, do.plot = TRUE, verbose.plot = FALSE)[["plot"]]
fit2_Dzp <- concRespCoreZR(row2_Dzp, do.plot = TRUE, verbose.plot = FALSE)[["plot"]]
fit3_Dzp <- concRespCoreZR(row3_Dzp, do.plot = TRUE, verbose.plot = FALSE)[["plot"]]
p1 <- cowplot::plot_grid(fit1_Dzp, fit2_Dzp, fit3_Dzp, nrow=2)
cowplot::plot_grid(SAplot_Dzp, p1, nrow=2)
Figure 6: BMC confidence interval overlap.
tcpl_out.dt <- as.data.table(DNT60_tcpl_out)
to.plot <- tcpl_out.dt[hitcall>0.8 & endp%in%c("avgS_L","avgA_L","avgJ_L","hbt_L",
"strtlA","strtlAavg","strtlF",
"avgS_D","avgA_D","avgJ_D","hbt_D",
"avgS_T","AUC_r")]
to.plot[, endp := factor(endp, levels=rev(c("avgS_L","avgA_L","avgJ_L","hbt_L",
"strtlA","strtlAavg","strtlF",
"avgS_D","avgA_D","avgJ_D","hbt_D",
"avgS_T","AUC_r")))]
to.plot[name=="Polybrominated diphenyl ether (PBDE)-47", name:="PBDE-47"]
ggplot(to.plot, aes(x=log10(bmd),y=endp)) +
geom_point() +
geom_errorbar(aes(xmin=log10(bmdl), xmax=log10(bmdu))) +
scale_x_continuous(breaks=c(-3,-2,-1,0,1,2)) +
facet_wrap(. ~ name)+
theme_bw() +
theme(text = element_text(size=14)) +
labs(title = "Overlap of BMC Estimates for Chemicals", x = "log(BMC)", y = "Endpoint")
Supplmental Figure 7: D-Sorbitol active curve-fits.
chemical <- "Amphetamine"
row1 <- rows_n[[chemical]][["avgS_L"]]
row2 <- rows_n[[chemical]][["hbt1_L"]]
row3 <- rows_n[[chemical]][["strtlF"]]
fit1 <- concRespCoreZR(row1, do.plot = TRUE, verbose.plot = FALSE)[["plot"]]
fit2 <- concRespCoreZR(row2, do.plot = TRUE, verbose.plot = FALSE)[["plot"]]
fit3 <- concRespCoreZR(row3, do.plot = TRUE, verbose.plot = FALSE)[["plot"]]
legend <- cowplot::get_legend(fit1)
title <- cowplot::ggdraw() +
cowplot::draw_label("Amphetamine Exposure: Active Curve-Fits",
x = 0, hjust = 0, size=24) +
theme(plot.margin = margin(0,0,0,7))
edit_plot <- function(plot) {
title <- plot$labels$title
plot$labels$title <- gsub(paste0(chemical," for "), "", title)
plot$labels$subtitle <- NULL
plot <- plot + theme(legend.position="none")
return(plot)
}
plots <- lapply(list(fit1, fit2, fit3), edit_plot)
p1 <- cowplot::plot_grid(plotlist=plots, nrow=2)
p2 <- cowplot::plot_grid(p1, legend, rel_widths = c(3, .4))
cowplot::plot_grid(title, p2, ncol=1, rel_heights = c(0.1,1))
Figure 7: BMC Heatmap
draw(htlist, merge_legend = FALSE, annotation_legend_list = lgd_list, annotation_legend_side = "top", align_annotation_legend = "heatmap_center")
Concordance Between this Study and Other Studies
How many chemicals had papers associated with them?
litReviewSummary[Number.of.associated.publications!=0, .N]## [1] 37How many papers of the 36 papers report chemical effect.
length( litReviewConc[neurotox==1 & DOI!="DNT60", unique(DOI)] )## [1] 30Supplmental Table 5: Number of Associated Chemicals and Reproducibility of Active Chemicals
litReviewSummary1 <- litReviewSummary[,-c(2,3,4)]
DT::datatable(litReviewSummary1, caption="Formatted Literature Review table.",
extensions = 'FixedColumns',
colnames = c("Chemical Name", "Number of Reviewed Publications", "Number of Reviewed Publications Reporting Effect", "Percent of Publications Reporting Effect"),
options = list(scrollX=TRUE),
rownames = FALSE)Find median LOELs for DNT60 study.
DNT60_medLOELs <- DNT60LOELs[, .(cpid, med_LOEL = matrixStats::rowMedians(as.matrix(.SD),na.rm=TRUE)), .SDcols=names(DNT60LOELs)[-1]]Table 6: Use median gabi LOEL to compare to other studies.
compareLOELs <-litReviewLOELs[DOI != "DNT60"][DNT60_medLOELs, on = .(cpid)]
setnames(compareLOELs, "med_LOEL", "DNT60_med_LOEL")
compareLOELs[, LOEL_diff := neuro_LOEL - DNT60_med_LOEL]
compareLOEL_summary <- compareLOELs[, .(mean_LOEL_diff = mean(LOEL_diff,na.rm=T), sd_diff_mag = sd(LOEL_diff,na.rm=TRUE)), by = .(cpid)][order(sd_diff_mag)]
compareLOEL_summary[, `:=` (mean_LOEL_diff = signif(mean_LOEL_diff,digits=3), sd_diff_mag = signif(sd_diff_mag,digits=3))]
DT::datatable(compareLOEL_summary,
caption="Summarize Differences in LOEL Magnitudes",
extensions = 'FixedColumns',
options = list(scrollX=TRUE),
rownames = FALSE)The Effect of Developmental Toxicity on the Detection of Neurotoxicity
How many chemicals have at least one associated paper that reports chemical activity?
litReview_actives <- litReviewConc[DOI!="DNT60" & neurotox==1, unique(cpid)]
length(litReview_actives)## [1] 31# Gather data for later plots and questions
DNT60LOELs.melt <- melt(DNT60LOELs, variable.name = "endpoint", value.name = "LOEL")
DNT60LOELs.melt[, `:=` (DOI="DNT60", category="intra_neuro")]
compareLOELs <- litReviewLOELs[DOI != "DNT60",-"dev_LOEL", with=FALSE][, .(cpid, DOI, LOEL=neuro_LOEL, category="inter_neuro")]
dev_LOELs <- litReviewLOELs[, .(LOEL=unique(dev_LOEL),category="intra_dev"), by=.(cpid)]
compareLOELs <- rbind(DNT60LOELs.melt, compareLOELs, dev_LOELs, fill = TRUE)gabi_LOEL_actives <- compareLOELs[category=="intra_neuro" & !is.na(LOEL), unique(cpid)]
discordant_inactives <- litReviewSummary[!cpid%in%gabi_LOEL_actives, .(cpid, Number.of.associated.publications, Number.of.publications.with.activity, `%.of.papers.with.activity`)][`%.of.papers.with.activity`>50][order(-`%.of.papers.with.activity`)]How many chemicals with a gabi behavioral LOEL had an associated publication report observed effect at at least one concentration.
table(gabi_LOEL_actives %in% litReview_actives)##
## FALSE TRUE
## 6 15What chemicals were active in other publications (>50%) but inactive in our analysis?
discordant_inactives## cpid Number.of.associated.publications
## 1: Bis(tributyltin) oxide 2
## 2: Cadmium chloride 1
## 3: Caffeine 1
## 4: Dieldrin 2
## 5: Isoniazid 1
## 6: Maneb 1
## 7: Nicotine 3
## 8: Phenol 1
## 9: Lead acetate 5
## 10: Acrylamide 3
## 11: Colchicine 3
## Number.of.publications.with.activity %.of.papers.with.activity
## 1: 2 100.00000
## 2: 1 100.00000
## 3: 1 100.00000
## 4: 2 100.00000
## 5: 1 100.00000
## 6: 1 100.00000
## 7: 3 100.00000
## 8: 1 100.00000
## 9: 4 80.00000
## 10: 2 66.66667
## 11: 2 66.66667discordant_inactives[Number.of.associated.publications > 1]## cpid Number.of.associated.publications
## 1: Bis(tributyltin) oxide 2
## 2: Dieldrin 2
## 3: Nicotine 3
## 4: Lead acetate 5
## 5: Acrylamide 3
## 6: Colchicine 3
## Number.of.publications.with.activity %.of.papers.with.activity
## 1: 2 100.00000
## 2: 2 100.00000
## 3: 3 100.00000
## 4: 4 80.00000
## 5: 2 66.66667
## 6: 2 66.66667Table 7: Median differences between developmental LOELs and study behavioral LOELs.
LOEL_diff_DN <- dev_LOELs[,.(cpid,dev_LOEL=LOEL)][compareLOELs[category == "inter_neuro", .(cpid, DOI, LOEL)], on = .(cpid)][, .(cpid, DOI, mag_diff=LOEL-dev_LOEL)]devBhvLOEL.diff <- LOEL_diff_DN[, .(median_mag_diff = median(mag_diff,na.rm=T), sd_mag_diff = sd(mag_diff,na.rm=T)), by = .(cpid)][order(-median_mag_diff)]
DT::datatable(devBhvLOEL.diff[!is.na(median_mag_diff), .(cpid = cpid, median_mag_diff = signif(median_mag_diff,digits=3), sd_mag_diff = signif(sd_mag_diff,digits=3))], caption="Median Difference in Magnitude from Developmental LOEL to Literature Review Behavioral LOELs",
colnames = c("Chemical Name", "Median Magnitude Difference", "SD in Magnitude Difference"),
extensions = 'FixedColumns',
options = list(scrollX=TRUE),
rownames = FALSE)Exclude those without majority of publications reporting activity.
cpid_maj_active <- litReviewSummary[`%.of.papers.with.activity`>50 & Number.of.associated.publications>1, cpid]
devBhvLOEL.diff[cpid %in% cpid_maj_active]## cpid median_mag_diff sd_mag_diff
## 1: Dieldrin 0.6365006 0.3373757
## 2: Nicotine 0.3124595 0.1118929
## 3: Heptachlor -0.1020600 0.7071068
## 4: Bis(tributyltin) oxide -0.2665894 1.0768379
## 5: Lead acetate -0.4586073 1.4744919
## 6: Colchicine -0.4900526 1.0098438
## 7: Valproate -0.6629268 0.3797085
## 8: Haloperidol -1.0401342 0.1938097
## 9: Bisphenol A (BPA) -1.4259687 0.5599005
## 10: Diazepam -1.6020600 0.0000000
## 11: Chlorpyrifos -1.6243364 1.1082443
## 12: Fluoxetine -1.9705933 0.5671248
## 13: Acetaminophen NA NA
## 14: Acrylamide NA NAFigure 8: Developmental and nuero- toxic LOELs observed in this study and others.
compareLOELs1 <- data.table::copy(compareLOELs)
cpid_wO_pub <- litReviewSummary[Number.of.associated.publications==0, unique(cpid)]
cpid_w_pub <- litReviewSummary[Number.of.associated.publications!=0, unique(cpid)]
# Compare LOELs has different spellings
# Bis(tributyltin) oxide
# Diethylene glycol
compareLOELs1[cpid%in%cpid_w_pub, facet := 1]
compareLOELs1[cpid%in%cpid_wO_pub, facet := 2]
# Order chemicals in each facet first by developmental toxicity and then by median LOEL were devTox doesn't exist and then by median study LOEL where median LOEL doesn't exist
factor.orderA1 <- rev( compareLOELs[cpid%in%cpid_w_pub & category=="intra_dev" & !is.na(LOEL)][order(-LOEL), unique(cpid)] )
factor.orderB1 <- rev( compareLOELs[cpid%in%cpid_wO_pub & category=="intra_dev" & !is.na(LOEL)][order(-LOEL), unique(cpid)] )
factor.orderA2 <- DNT60_medLOELs[!(cpid%in%factor.orderA1) & cpid%in%cpid_w_pub & !is.na(med_LOEL)][order(med_LOEL), unique(cpid)]
factor.orderB2 <- DNT60_medLOELs[!(cpid%in%factor.orderB1) & cpid%in%cpid_wO_pub & !is.na(med_LOEL)][order(med_LOEL), unique(cpid)]
factor.orderA3 <- compareLOELs[!(cpid%in%factor.orderA2) & !(cpid%in%factor.orderA1) & cpid%in%cpid_w_pub & category=="inter_neuro", .(med_LOEL = median(LOEL,na.rm=TRUE)), by = .(cpid)][!is.na(med_LOEL)][order(med_LOEL), unique(cpid)]
factor.orderB3 <- compareLOELs[!(cpid%in%factor.orderB2) & !(cpid%in%factor.orderB1) & cpid%in%cpid_wO_pub][order(-cpid), unique(cpid)]
factor.orderA4 <- compareLOELs[!(cpid%in%factor.orderA3) & !(cpid%in%factor.orderA2) & !(cpid%in%factor.orderA1) & cpid%in%cpid_w_pub][order(-cpid), unique(cpid)]
factor.order <- c(factor.orderA4, factor.orderA3, factor.orderA2, factor.orderA1, factor.orderB3, factor.orderB2, factor.orderB1)
compareLOELs1$cpid <- factor(compareLOELs$cpid, levels=factor.order)facet_labels = list("1" = "Chemicals with Associated Publications",
"2" = "Chemicals without Associated Publications")
facet_labeller <- function(variable, value) {
return(facet_labels[value])
}ggplot() +
geom_point(compareLOELs1, mapping = aes(x=LOEL, y=cpid, color=category, shape=category), size = 3) +
scale_shape_discrete(labels = c("Extra-Study Behavioral LOELs","Developmental LOEL","Study Behavioral LOELs")) +
scale_color_manual(labels = c("Extra-Study Behavioral LOELs","Developmental LOEL","Study Behavioral LOELs"), values=colors[c(1,8,4)]) +
theme_bw() +
theme(text = element_text(size=14)) +
facet_wrap(facet~., dir = "v", as.table = TRUE, drop = TRUE, scale = "free_y", strip.position = "right", labeller = facet_labeller) +
labs(shape = "Critical Doses",
color = "Critical Doses",
x = paste0("Concentration (log10(",paste0("\U03BC","M"),"))"),
y = "Chemical Name")
Supplemtnal Figure 8: Unity plot of median potencies.
# plot medians with quantiles and interquartile range
cHairPlot.data1 <- compareLOELs[category=="intra_neuro", .(med=summary(LOEL)[3], Q1=summary(LOEL)[1], Q3=summary(LOEL)[5]), by=.(cpid)]
cHairPlot.data2 <- compareLOELs[category=="inter_neuro", .(med_x=summary(LOEL)[3], Q1_x=summary(LOEL)[1], Q3_x=summary(LOEL)[5]), by=.(cpid)]
# _x indicates values are derived from "extra" laboratory studies
cHairPlot.data <- merge(cHairPlot.data1, cHairPlot.data2, all = TRUE)ggplot(cHairPlot.data, aes(x=med, y=med_x)) +
geom_point() +
geom_abline(intercept=c(-1,0,1), slope = c(1,1,1), linetype=c("dashed","solid","dashed")) +
ggrepel::geom_text_repel(aes(label=cpid,color=cpid), box.padding = 0.5) +
lims(x=c(-2,2), y=c(-2,2)) +
labs(x = "Median Neurotoxic LOEL Observed in this Study",
y = "Median Neurotoxic LOEL Derived from Literature Review Studies",
title = "Concordance of Neurotoxic LOELs Across Studies") +
theme_bw() +
theme(legend.position="none")